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Abstract 

Various approaches to the numerical representation of the Incomplete Gamma Function F m (z) 
for complex arguments z and small integer indexes m are compared with respect to numerical fitness 
(accuracy and speed). We consider power series, Laurent series, Gautschi's approximation to the 
Faddeeva function, classical numerical methods of treating the standard integral representation, 
and others not yet covered by the literature. 

The most suitable scheme is the construction of Taylor expansions around nodes of a regular, 
fixed grid in the z-plane, which stores a static matrix of higher derivatives. This is the obvious 
extension to a procedure often in use for real- valued z. 
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I. OVERVIEW 



Motivation 



The Incomplete Gamma Function F m (z) is at the heart of the computation of Electron 
Repulsion Integrals over Gaussian-type basis functions l[ 3, If these are attached to 

moving atoms, and the dominant part of the time dependency is kept with the bases (instead 
of being hidden in the expansion coefficients), the argument z is complex- valued (App. |A"|) . 

With respect to the index m, we look at this as being derived from coupling of integer- 
valued orbital quantum numbers, whence deal only with small, non-negative m unless oth- 
erwise noted. 

Several generalizations □ □ □ are not considered here, except the — important- 
identification with some confluent hypergeometric series. 



B. Contents 



The explicit intention of this script is to compare a wider range of methods than proposed 
on the same subject before M. Continued fractions are not covered here, because they have 
a, ready been detaUed beforelQ. Mso, th e Te m ,ne and Pan, app—ions flS have 
been left aside, because they would work around the Complementary Error Function of a 



complex argument, the calculation of which is already of the same order of complexity as the 
original problem. No attempt is made to discuss the Gauss-Rys quadrature jl3 
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for complex z, which constructs a system of orthogonal polynomials over a finite interval 
with weight function exp(—zt 2 ); the explicit notation of its polynomial of degree 1, Eq. 
(2.4) by Steen et al. 17| . illustrates that this also would start from the complex Error 
Function. Barakat [la} 19, (13.3.9)] reports on z-values on the imaginary axis mapped on 
Bessel functions with real-valued arguments, which we do not follow on the same reasoning. 
For the case of real-valued z and large m we refer the reader to the article by Takenaga 

□ 

F m (z) is eventually defined by an integral with a rather simple kernel. The following 
chapters are roughly grouped according to how much effort is spent on isolating special 
aspects of this kernel in exactly integrable terms, in hope of catching up the oscillations 
induced by Qz or the steep slopes induced for large 9ftz, for example. 



C. Fundamentals 



We define F m (z) through its simplest integral representation, 

F m (z) = f t 2rn e~ zt ' 2 dt = - f u m - 1/2 e- zu du. (1) 
Jo 2 J 

Though the application in the literature often focuses on the positive real axis (as z = 
x + iy represents some square of real- valued widths of orbital exponentials then), we will 
cover the general case. 

Anyway, the complex-conjugate symmetry 



F m {z) = F m {z), (2) 
(with z ee — iQz) allows us to restrict the analysis to the cases of Qz > 0. 



The forward recurrence 



2l|[19l (6.5.21)] 



2zF m (z) = (2m - l)F m ^{z) - e~ z (m > 1) (3) 

is a rephrasing of |3, (13.4.4)]. It suffers from cancellation of about — log 10 \2z\ (m = 0) 
and — log 10 \z/m\ (m > 0) decimal places when applied to \z\ 1 22|. The corresponding 
backward recursion suffers from cancellation of about — log 10 \ (m — l/2)/z\ decimal places, if 
$tz is negative and \z\ large. Generally speaking, the recursion allows to pick a numerically 
favorable m and switch to others at a small additional cost. 

Further down, numerical precision is demonstrated by the number of decimal places, 
defined as the negative Brigg Logarithm of the relative error, d = — log 10 |1 — F m (z)/F m (z)\, 
given a high precision actual value F m (z) and its approximation F m (z), both computed with 
Maple at 30 decimal places. 



II. HYPERGEOMETRIC SERIES 

Expansion of the exponential in ([T} and term-by-term integration yields the confluent 
hypergeometric series (Kummer Function) 

Fm(z) = 2mTT lFl(m+ ^ m+ | -2;); (m> " 1/2) (4) 

X F X (a-a + l-z ) = 1 + __ + ... + ___ + ... 5 

a + 1 (a + n)n\ 
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We start with the hypergeometric power series as a reference because there is an evident 
and flexible implementation with a clear numerical cost of one complex multiplication (which 
could be implemented with three real multiplications and five real-valued additions, or with 



four real-valued multiplications and two real-valued additions [23J, 3.7.2]) and addition per 
term: one accumulates terms until the new term's contribution falls below a limit set by a 
preset relative error. Fig. ^ verifies that the series converges best close to the origin of the 
complex plane, as expected for any power series. The unexpected feature is that the series 
performs worse if terms are alternating — that is if z is close to the positive real axis in Fig. 
[0 — than in the case of the non-alternating mirror point —~z. This is a by-product of a growth 
of the terms modulus up to the partial sum of index n « \z\ — 1. The alternating case must 
overcome a massive cancellation of digits when passing this index. Consequently it needs 
much more terms until the partial sums approach the order of magnitude of the exact result. 
If measured in terms of the absolute error after summation of n terms, the alternating case 
would indeed perform better. 

The Shanks transformation j^J €2(S n ) of the partial sums S n of © — which is to this 
second order just the Aitken transformation [lj], (3.9.7)] — would improve the accuracy of 
the plots of Fig. by roughly 1.5 digits. This involves handling of finite differences between 
numbers that are (supposedly) already close to each other and is more tricky than the 
analytical transformation formulas indicate. 

III. LAURENT SERIES 
A. Barnes' Analysis 

The asymptotically convergent Laurent series for large \z\ is a special case of Eq. (13.5.1) 
in ^| or taken from §6 in j^j: 

ae z 



iFi(o; a + l;z) = T(a + l)(-z)- a + — ^(1 - a) n z' n , (6) 
with Pochhammer's Symbol defined as 19, (6.1.22)] 



z 

n=0 



1 , (n = 0) 

(b) n ^T(b + n)/T(b)={ A ] (7) 

b(b + 1)(6 + 2) • • • [b + n - 1) , (n > 0) 
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FIG. 1: Contour levels of the number of valid decimal digits d = 4. . . 18 of F m (z) by the power 
series (@J and ©, if the power series is truncated after n = 30 or 60 terms. 
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FIG. 2: Contour levels of the number of valid decimal digits g? = 2 ... 12 of F m (z), if the Laurent 
scries © is summed up to n = \z\ + a. 



221. 



This series is also known under the label "high-T" expansion in quantum chemistry |4 , 
The Gamma Function is not of concern since it is only needed for half-integer values, and 
would be tabulated based on (6.1.12)]. Asymptotic convergence means that the terms in 
© shrink until n < \z\ + a and grow afterwards. This inherent limitation to the achievable 
accuracy is put into concrete with Fig. El 

Swapping the the sign of z in Eq. (jlj), the series becomes alternating near the positive 
real axes of the plots in Fig. |2 for 9ftz > 0, which leads to some obvious left-right asymmetry 
in the precision attained. 
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The contrasting regions of good convergence manifested in Fig. El and Fig. ^ suggest to 
combine these results into Fig. El The maximum number of terms needed this way to obtain 
d = 12 digits for m = is n < 92 for the entire z plane (Fig. El top), to obtain d = 12 digits 
for m = 1 is n < 84 for the entire z plane (Fig. El middle), and to obtain d — 17 digits for 
m = is n < 166 for the entire z plane (Fig. 01 bottom). 
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FIG. 3: Contour levels of the number n of terms needed to obtain an accuracy of d = 12 or 17 

digits of F m (z), if the power series © and the Laurent series © are used compliment arily. 
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B. Gargantini-Henrici Converging Factors 



The Gargantini-Henrici analysis 26[ of the converging factor of the Laurent series © 



allows a more accurate calculation of the truncated series 

oo 

'^Esr ( 8 ) 



^ oo oo 

E(l - a) n z~" X 

n=0 n=0 



(n) 


_ c n+l 


1 - 


a + n; 


(n = 0,l,2,...) 


( n 






(n) 


= k; 


[n = 


0,1,2,. 


■■) 


In) 


= n + k - 


- a; 


(n = 


0,1,2,...) 



The coefficients c n = (1 — aL are fed into the quotient- difference scheme easily derived from 
Sec. 5] or taken from 23, (3.9.3)], 

(9) 
(10) 
(11) 

and (JBJ) is approximated by 

AT-l 
n=0 

z a {N) p {N) a (N) p (N) 
# N {z) = JL?l-ZL-«Z-h-... (iV = 0,l,2,...) (13) 

The following results of Fig. H] are based on a "best knowledge" approach in the sense that 
the approximation (fT2"|) sums to the same N < \z\ + a as in the previous section, and that 
the continued fraction are accumulated until the lower indexes in and ejj. have reached 
iV — so to recycle the same c n that appear in the main series. (This roughly triples the 
number of multiplications and additions for a particular z compared to the approach of just 
truncating ©.) 

Compared with Fig. the multiplication with the convergent factor has approximately 
doubled the number of valid digits in a range of intermediate \z\, but as the rational function 
introduced by the continued fraction has been allowed to grow to polynomial degrees of 
numerator and denominator comparable to the cut-off order of the series, there is no longer 
a monotonic increase of accuracy away from the origin. The regions in the complex plane 
of predictable accuracy have got a complicated shape. 
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This representation gets slightly more compact by decomposition into partial fractions and 
use of the symmetry n <-> n — k in the fc-sum: 

F 2 fs = t-*J a / 15) 

ml ; ^2m + n + l^A;!(n-A;)!(2m + 2A; + l) 1 ; 

n=0 fc=0 V ; V y 

The decrease of the coefficients of this power series is demonstrated in Tab. UJ The 
competitive power series (0) converges faster since its coefficients decrease faster, roughly 
oc l/(n • n!) as a function of n. 



B. Power Series with Split-Off Exponential 

The fundamental power series (jlj is derived from replacing e~ 2 * 2 by its power series; by 
any truncation of the series, this represents the integral kernel just at t = 0. We investi- 
gate a more accurate interpolation, which equals the kernel at both limits of the t-interval 
[0, 1], which splits off a simpler exponential that still can be integrated exactly, and which 
accumulates the (smaller) remainder in a different power series: 

e -zP = e -zt + z[t _ e) _ Z\ t 2 _ t 4 } + £_\ t 3 _ t 6) _ . . . (16) 

Fm(z) = / t2me ' Ztdt ~ E (i-l)!(2m + n + l)(2m + 2n + l)- 



In particular, 



^ (n - 1)! (2m + n + l)(2m + 2n + V 

Fo(z) = ~z~ ~ ^ (n - 1)! (n + l)(2n + 1) ' (18) 
As a side note, a decomposition in partial fractions and insertion of the hypergeometric 
notation with (@J) at m = 1 yields 

ov ' z ^ n! n + 2)(2n + 3) v 7 

n=0 v y v ' 

^""l.-f ( "'''" 1 -- V^— (20) 



2! ' n\ n + 3/2 ' n\ n + 2 

n=0 ' n=0 



§l.Fi(3/2;5/2;-2)=2Fi(2) | 1 F 1 (2;3;-z)=[l-(z+l)e-^ ]/2 2 

and emerges as a complicated adornment of 0. Supposed one has a fast, reliable method 
to compute 1 — e~ z , (fT7j) looks beneficial compared to (j5J) because the total power of the 
summation variable n in the denominator is slightly larger. Graphing the results of (j!8)) the 
same way as in the two upper plots of Fig. ^ however, would yield no differences visible to 
the eye. Therefore we do not look into this ansatz further. 
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16 (-10)0.11574 62528 45873 78175 30636 17 


(-11)0.73931 27958 11477 67351 64672 14 


17 (-11)0.12095 88230 02785 23102 13378 56 


(-12)0.79139 24236 39157 01672 60238 06 


18 (-12)0.12018 49078 18966 37363 41310 02 


(-13)0.80345 79062 12657 33014 90722 55 


19 (-13)0.11380 77808 88498 70226 59025 07 (-14)0.77576 42229 41899 25361 71656 84 


20 (-14)0.10293 02165 82138 49790 10396 95 (-15)0.71410 21050 93626 04408 94055 80 



TABLE I: Decrease of the coefficients in ()15p . cases m = and m = 2, as a function of n. 
The numbers in parentheses denote multiplication by powers of 10, e.g., for n = 17, we have 
0.1209 . . . • 10~ n and 0.791 . . . • 10~ 12 . 

C. Power Series of the Half Argument 

Convergence of power series is generally faster closer to the origin; the trigonometric 
identity sin 2 (u;/2) = [1 — cosu;]/2 allows us to reduce the distance between z and the origin 
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by half if we substitute t = sm(u/2) in (JTJ): 



Fo(z) = - I e- zt dt 



-z/2 PIT 



OJ 



e 2 COSUJ cos-duj 



e z l 2 ^ 1 fz\ n 



J ^ n\ V2 

n=0 



n W J 

cos oj cos — au;. 



The auxiliary cu-integrals would be drawn from the recursion [23, (2.538.1)] 

1 



OJ 



cos n cos — dw = 

2 n + 1/2 



n / cos™ 1 cos — gL> 



(21) 



(22) 



Since these are of the order of 1 for all n, we are left with a power series which converges 
~ z n /(n\2 n ), which is to be compared to ~ z n /(n\n) of (J3J). 
The generalization to nonzero m reads 



F m {z) 



e'*/ 2 ^ 1 fz 



Y- - 

4 ^ n\ V2 

n=0 



sin m — cos" a; cos — doj 



(23) 



The (n, m)-table of the tu-integrals could be generated from the table at m = or from 
scratch, 



• 2m ^ n 

sin — cos oj cos 
2 



—doj = — Y(-) k ( m I / cos fe+n cos -oL> 
2 2 m ^ v J \k J J n 2 



fc=0 
m n 

EE 



m \ I n 



ik+n+lnl+l (^ + 0' 



(24) 



fc=0 z=o 

The examples of Fig. 03 demonstrate that about 8 additional digits have been gained for the 
case m = and n = 30 terms relative to the data of Fig. ^ There is no additional run-time 
cost since the coefficients table of (|2*3j) is static without z-dependence. 
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FIG. 5: Contour levels of the number of valid decimal digits d = 4. . . 20 of F m (z) by the power 
series (|'23|) . if the power series is truncated after n = 30 or 40 terms. For m = and n = 60, d is 
> 19.6 in this z-domain. ' ' 



V. INTERPOLATING THE INDEX 



The integral is solvable if 2m = 2a — 1 is an odd integer |28|, (2.321.2)]: 

,(a = 1,2,3,...). (25) 

This could also be derived from (j^J) and (|fi|). where the sum in ((HJ) terminates if a is a positive 
integer. 

Given z, let (|23j) be computed for an index set m — 1/2, 3/2, 5/2, . . . ,N+ 1/2, optionally 
facilitated by the index recursion j3, (2.321.1)] or (jSJ). Let us pursue the idea that approx- 
imate values at intermediate m — 1, 2, . . ., which we are actually interested in, are deduced 
by some interpolation. 

If this is done by the unique interpolating polynomial J^l bjvn? of degree N in m, the 
intermediate step of the calculation can formally be written down as a (N + 1) x (N + 1) 
inhomogeneous system of linear equations to get the N + 1 unknown polynomial coefficients 
bj (a Lagrange Interpolation might be cheaper numerically, though), 

TV 

J2Kbj=F mk (z), (k = 0,1, ... ,N;m k = k + 1/2), (26) 

j=0 

and we get for example Fig. El We see: (i) The method becomes more and more unreliable 
as Jfc increases, which emerges from different weighting by the factors t 2m and e~ zt2 in the 
integral kernel: The derivative with respect to m multiplies the integrand (0) by 2 Int. To 
keep the derivative small relative to F m (z) itself (and to keep F m (z) a flat function of m), 
it is advantageous that t 2m e~ zt2 weights stronger at the right limit of the t-interval, where 
hit stays small. t 2m becomes large near t = 1, whereas e~ zt2 is larger near t — or t — 1 
depending on the sign of 3tz.) (ii) The estimates are better close to the middle of the m 
interval [|, iV+1], which was sampled to define the interpolation polynomial, than for m close 
to the interval limits — which is expected for any interpolation derived from approximately 
equidistant sampling points. 

No such interpolation polynomial of m satisfies (JBJ) on the entire real m-axis — though 
the exact F m (z) does for all m > 1/2. Therefore one could seek after improvement of this 
interpolation by enforcing some compliance with (jSJ). Test calculations were made following 
the strategy that one or more lines in the system of linear equations ()26|) are replaced by 
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coupling some of the unknown F m (z) via (j3J), in concrete 

N 

J2[^mi-{2m u -l)mi}b j = -e- z , (m v = 0, 1, 2, . . . ; m u = m v + 1). (27) 

3=0 

The polynomial that ensues does no longer hit the F mk (z) that was thrown out, but follows 
(jSJ) for those m u brought in. Results of this ansatz with N = 4, removal of Fg/2 from the 
set of interpolated points, and introduction of the m v — 1 m u = 2 coupling with (|2~7jl 
add at most about half a digit of accuracy to what is already shown in Fig. |U1 There is no 
further improvement if the line for F7/2 in (|27)|) is also removed to add the line (|27j) for the 
m v = 2 <-» m u = 3 coupling. 

Furthermore one would try to add known information on the first derivatives with respect 
to m to enhance the quality of the polynomial interpolation. Unfortunately, even the simplest 
dF m (z) / dm at m = 1/2 would demand computation of the Exponential Integral 

dJ ^- =~ [E l {z)+ 1 + \og{z)\. (28) 
dm \m=i/2 2z 

This infects also the partial derivatives at m — 3/2, 5/2, . . . 

^F m (z) = -F m ^(z) + —-^—F m _i(z); (m > 1/2). (29) 

am z z 8m 

Another attempt of refinement is to acknowledge the simple poles at the negative 
half integers of m. Test calculations with the modified separation ansatz 30( F m (z) « 
(Sjlo bjmi) I (2m + 1), which manifests the pole at m = —1/2, and again a numerator poly- 
nomial of degree N = 4 result in changes of up to one digit (in both directions) compared 
to the bottom graph in Fig. 

In summary, it seems to be difficult to bridge the gap between F m at half integer and 
integer m through interpolation. 



case o 
Ei(z) 
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VI. GENERIC METHODS OF INTEGRATION 



A. Local Taylor Expansions in the Integration Interval 

1. Expansion of the Exponential 

The "global" first approximation of the integrand in Sec. UVBl mav be pushed one notch 
further towards a brute-force numerical method by slicing the t-interval [0, 1] into N same 
size subintervals of half-width A = 1/(2JV), centered at t, = (21 - 1)/A (I = 1, . . . , N). In 
each of these subintervals, exp(— zt 2 ) is approximated by its Taylor series around t\. The 
derivatives are [3, (7.1.19)] 

{jX^ = (-) n z n/2 H n (^-zt)e- zt \ (30) 

in terms of Hermite Polynomials H n — looking at the exponential as the first derivative of 
the error function — , whence the Taylor series 

oo 1 

e- zt ' 2 =J2^(t- Un-rz^H^t^e-^. (31) 

n=0 

Fq(z) is the Riemann sum over the subintervals 

N pti+A N 00 An+1 

f (z)^J2 dte ~ zt2 = 2 E e ~"* E rrw /% »(^) ( 32 ) 

l=l Jt l- A 1=1 n=0,2,4,... [U+ j ' 

Fig- HI shows the accuracy of ()32jl for two AT, keeping the polynomial expansion only up to 
some degree n. 

The two graphs in Fig. [7|with degrees kept up to n = 6 show the two spots at z ~ 2.80 
and z ~ 7.02 with higher precision than their surroundings. This is part of a more general 
phenomenon, in which for some z the positive and negative lobes of the Hermite Polynomials 
are sampled with best effective cancellation. This generates one such spot at z w 4.08 
(n < 4), those two for n < 6,. . . 

On the computational expense: as the formula requires Hermite polynomials of even 
indexes only, there is no need to compute the yfz and no gain using their recursion formulas. 
Two (complex) multiplications compute zt 2 . Each H n is a polynomial of degree n/2 in 
this combined variable, which costs n/2 multiplications and n/2 additions with the Horner 
scheme. The powers of A are fixed, and there are about 3 multiplications for each term in 
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FIG. 7: Contour levels of the number of valid decimal digits d = 8 . . . 16 of Fq(z), if the series (|32|) 
is used with N = 10 or 20 subintervals and if the Hermite Polynomials are included up to degree 
n = 6 or 8. Note that in a region of Jftz <J 8 aric? small |3f^| an increase of n from 6 to 8 decreases 
the acciiracv. 



n, plus one (complex) exponential. This is to be multiplied by N. This totals at least 16iV 
multiplications if n runs up to 6, and so only the lowermost picture in Fig. [7| would be part 
of a fair comparison with the other competitive approaches. 



2. Expansion of the Algebraic Factor 

In a similar manner as above, one could expand the algebraic factor u m ~ 1 ^ 2 of (JIJ in 
Taylor series around u\ = (21 — 1)/A to end up with closed form integrals. To second order, 

m-l/2 _ m-1/2 . / 1\ m-3/2/ s . / 1\/ 3 m-5/2 (u ~ Ui) 2 

u 1 +(m--)u l ' (u-uj) + (m--)(m--)u, , (33) 

and the Riemann sum 

F m (2) w -— V e-^ M r 1/2 (e- zA - e zA + m - 1/2 [ e - zA (l + zA) - e zA (l - 2A 

1=1 1 

(m- l/2)(m-3/2) r ~ 2a2 ~ 2a2 



+ 



A + ^ A + l)_ e .A ( l^__ JgA+1 J 



v 2 ' v 2 

This achieves up to d = 2.6 decimals (m = 0, iV = 20), d = 5.4 (m = 1, N = 20), d = 5.8 
(m = 1, N = 40), and d = 8.2 (m = 2, iV = 40) in the z-domain as in Fig. (These 
numbers refer to the "lower left" corner of the z-region, and are a few digits worse in the 
opposite corner.) 

The main obstacle to higher performance is the poor fit of ()33|) close to u = 0. We may 
patch this by replacing the contribution in this subinterval, the term I = 1 where < u < 2A 
remains small, by the associated power series of the exponential, 

1=2 n=0,l,2,... V 1 / 1 / 

With this ansatz and the sum over n kept up to n = 3, the maximum number of digits in 
the z-domain as in Fig. [7|rise to d = 7.8 (m = 1, N = 40), and d = 10.4 (m = 2, N = 40). 

B. Fourier Expansion of the Algebraic Integral Kernel 

A Fourier expansion of the algebraic term in (JIJ 

00 

^qcos(Zw-) (35) 



u m-l/2 ^ ■ " 



1=0 
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m N 

512 256 128 

1 (-2)3.6 (-2)5.1 (-2)7.4 

2 (-5)4.6 (-4)1.3 (-4)3.7 

3 (-6)8.2 (-5)3.5 (-4)1.4 

TABLE II: Typical maximum deviation between the left hand side and the right hand side of ()35|) 
in the interval u £ [0, 1], if one period of f(u) is sampled by N points and the sum (|35|) limited to 
< / < N/2. Numbers in parentheses are powers of 10 as in Tab. EJ 



offers the series 



where c = 
oscillations 



and c 2 = c 4 = c 6 = . . . = have already been assumed. To reduce any Gibbs 



3jj of (|33|) at the ends of the interval [0, nto the even, 

4-periodic, and steady carrier function f{u) = u 171 ^ 1 / 2 for u E [0, 1], f(u) = 2 — (2 — n) m_1 / 2 
for u G [1,2], f{u) = f{2 — u) and f{u) = f(u + 4) elsewhere. (Obviously, the singularity 
at u = reduces the quality of this approach right from the start if m = 0.) The q are 
approximated by a discrete cosine transform on iV grid points 

N/2-1 



Cj N 



^ f(4k/N) cos(27r^) 

k=l 



(j = 1,2,..., N/2), (37) 



and the summation (j3T?j) is truncated at I = N/2. The q would be kept in constant tables 
since they do not depend on z, and the cost of evaluating (|3l)|) amounts to about N/4 
evaluations of the rational term. Fig. |S] shows that for moderately small m a precision of 
just of the order of 5 digits result from this type of evaluation, which is not efficient compared 
to other methods proposed here This is ultimately a progression of residual fitting errors 
(Gibbs oscillations) which remain rather large close to u = and u = 1 (Tab. ILT|) . 

In the special case where one would like to tabulate F m (z) along some fixed ordinate 
5fc = const in equidistant steps of Qz, one could draw a lot of additional benefit from the 
inherent parallelism of Fast Fourier Transform techniques applied to the decomposition 

1 f 1 

F m (z) = - / u m - 1/2 e- Uzu e- &zu du. (38) 

2 Jo 

This deems to be too special to be put into detail here. 
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FIG. 8: Contour levels of the number of valid decimal digits d = 1.5 ... 8 of F m (z), if the Fourier 

series (|Mfi|) is used with N = 256 or 512 abscissa points of the Fourier interpolation. 
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C. Taylor series in the complex z-Plane 

With [3 (13.4.9)], the complex derivatives of F m (z) are equivalent to ladder-type oper- 
ation with respect to m and a: 



n 



— ) iFi(a; a + 1; z) = — - — iiq(a + n; a + 1 + n; z), (39) 
dz J a + n 

" F m (z) = (-) n F m+n (z) (40) 



This close link between the index and the higher derivatives means that one may tabulate 
the expansion coefficients of the Taylor series 

e 2 (-e) n 
F m (z + e) = F m (z ) - eF m+1 (z ) + ^F m+2 (z ) + ■■■ + ^-F m+n (z ) + ■■■ (41) 

anchored at some z in the complex plane for all m at the same time. This keeps these 
tables smaller than for any function with decoupled index and argument. Despite the fact 
that these tables need to contain complex values and need to be arranged on a 2D grid in 
the complex plane, there is no inherently new aspect over what is already assessed in the 

n 

literature for real z [21] . Chebychev approximations in the complex case could be derived by 
expanding (— e) n in a sum over products of 5te and normalization to the interval [—1, 1] 
and independent progression for the real and imaginary part as described in [19I §22.20]. 

Fig. El is an example where the nodes span the Jizo interval from —33 to 18 with a stride 
of s = 3, and ^Zq the interval from to 36 also with a stride of s = 3: 

z = ks + ils (ife = -11, -10,..., 6; J = 0,1,..., 12) (42) 

Up to 23 terms need to be accumulated in (J41|) to calculate F (z) if z falls inside this finite 
grid's domain; values outside are handled with Eq. |H1 This term count is stored as entry 23 
at m = and d = 14 in the upper part of Tab. IIII1 A free entry in the table indicates, what 
could not be computed if the tabulated Fj(z) would be limited to j < 30, which would be 
less than m + n for this slot. The lower part of the table illustrates, how a three times denser 
grid of nodes cuts down on the worst case convergence, as it reduces the maximum distance 
|e| to the nearest zq from 3/a/2 to l/y/2. This necessitates a ninefold larger static table, 
unless one uses the smaller number of terms to reduce the maximum index j of tabulated 
F r 
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n 




FIG. 9: The number n of terms in the Taylor series Eq. (|41|) to achieve d = 14 digits of accuracy 
of F (z). 

m d 

12 14 15 16 17 

21 23 24 25 26 

1 21 23 24 25 26 
3 21 23 24 25 26 
5 21 23 24 25 

m d 

12 14 15 16 17 

14 16 17 17 18 

1 14 16 16 17 18 
5 14 15 16 17 18 

TABLE III: The upper table shows the maximum number of terms needed for nearest neighbor 
Taylor expansions (|41[) as a function of m and accuracy d on a support grid with spacing 3 in the 
z plane. The lower table illustrates the case of a denser spacing of 1. 
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On a side note, (j4"Tj) could be rewritten similar to |32j, (2.1)] 

F m {z Q + e) = e- Z0 F m (e) + z V — — 1 / F m+1+n (z ), (43) 

: ' m + 1/2 + n n! 



71=0 



using (J3J) once for each term, or integrating (JTJ) by parts. 
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D. Gauss-Jacobi quadrature 

As considered by Gautschi (3^] for a more general case, (JTJ) is readily accessible by a 
Gauss-Jacobi quadrature, 

n 

F m (z)~Y. w * e ~ Ztl - ( 44 ) 
i=i 

Weights Wi and abscissae tj are discussed in App. |EJ 

Fig. El demonstrates that this method is most robust at small \z\, which is expected since 
the Gauss quadrature effectively approximates exp(— zt 2 ) by a polynomial of degree 2n in t. 
The numerical expense roughly adds up to n computations of exponentials exp(— zt 2 ). This 
is cheaper than evaluation of (J32*j) at the same N = n; by further comparison of Fig. ITU1 with 
Fig. we conclude that the Gauss-Jacobi ansatz proposed here is superior to the method of 
Sect.lVDD 

The drawback by further comparison with the reference calculation of SeclHlis that each 
complex exponential needs much more CPU time than a complex multiplication. Tests with 
the C++ implementation by the Sun Forte Developer 7 Collection suggest a factor of about 
twenty. 

This is also the reason why the trapezoidal rule, and higher rules like the Simpson Rule 
which follow from a Richardson extrapolation, have been kept aside in this manuscript. 



E. Cubic Spline Interpolation 

The exponential of (JTJ) could be approximated by cubic splines in N subintervals [tj,tj + 
1/iV] (j — 0, . . . , N — 1) which cover the /-interval, to yield a sum over elementary integrals, 

F m ( z ) f2m h + Cit + C2f2 + Cst3 ] dt ( 45 ) 

j=0 A 

In each of these intervals, the four coefficients c« are defined by demanding that the cubic 
polynomial fit and its first derivative equal the exponential and its first derivative at both 
limits, tj and tj + %: 
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FIG. 10: Contour levels of the number of valid decimal digits d = 3 ... 19 of F m {z) by 1)44 jl . if the 

order of the Gauss-Jacobi quadrature is n = 10 or 20. 
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FIG. 11: Contour levels of the number of valid decimal digits d = 1 ... 8 of F m (z) by the spline 
interpolation (|45|) for N = 10 subintervals. One would gain one digit throughout this z region by 
increasing N to 20. The accuracy for m = 3 would be a few digits lower. 



The number of multiplications to solve this system of linear equations in each interval 
looks prohibitive, even though one would recycle the matrix elements and right hand sides, 
and even though the matrix is already close to triangular form. One must evaluate N 
exponentials and insert the q into about 2N polynomials of degree 2m + 4 to finalize (|45jh 
Actually, inserting (|46|) into ()45|) yields rather tight formulas (see j^] and App.[OJ), namely 

(47) 



Fo(z) 



j=0 



and 



AT-l 

E 



30 



(2zi 4 +1 - zt 2 J+l t 2 j + 8t 2 j+1 + bt j+1 tj - zt j+1 t) + 2tJ) e-^+i 



(4£ 



The drawback formulated in Sec. IVI Dl however, remains: we consider only very small N, 
and conclude from Fig. [n] that this numerical expenditure is too high to consider this method 
a competitive candidate. 
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F. Salzer's Numerical Inverse Laplace Transform 



F m (z) can be written in terms of the inverse Laplace Transform of some function P(m + 
1/2, z) which is loosely related to the v 2 probability distribution 

F ra (z) = 5p( fl ,z); (a = m + l/2), (49) 

where 

i pc+ioo i i 

p(a ' z) = 2m /_ (0O e "lJFTTr ds '' (c>0) ' (50) 

The binomial expansion of (s + l) a in powers of 1/s and interchange of integration and 
summation transforms (J5Uj) to the power series defined by (@J) and (JSJ). 

Substituting zs = p in P(a,z), ()50|) is in shape for Salzer's 35j approximation of the 
kernel l/[s(s + l) a ] by polynomials in 1/p. 

()49j] and fl5()l) assume the existence of the Laplace Transform of z a F m (z); this limits this 
proposal in general to ^kz > as indicated by the factor e z in (jBJ) and made explicit in [36, 
p. 341]. Therefore we use [11 1 

i rc+ioo i i 

P(a,z) = 1-Q(a,z) = 1 -ae~ z z a / e as — ds (51) 

to complement for $lz < 0. This, however, is of purely experimental nature since (|51|) 
depends on a contour integration passing between the pole and the branch point, and this 
is not at all accounted for in Salzer's sampling of the complex p-plane. 

Sample outputs of this approach are gathered in Fig. [12] using abscissa and weights as 
tabulated in Tab. IIV1 The discontinuity in the graphs at passing the imaginary z-axis is 
due to switching between (joTi|) and (joTj) . We see that the accuracy of Fi(z) rises by about 
one decimal if the number n of evaluations of the kernel is increased from 16 to 30, but the 
rise in the case F^{z) — not demonstrated in the figure — is only about a third decimal. The 
major difference to the Taylor and Lauren series methods of sections [H] IIII Al and IVI CI is 
that this here still demands computation of a complex valued root for each of the sampling 
points, and that prediction of the accuracy both as a function of m and z is complicated. 

Finally, we did not try to apply this approximation to the Laplace transform of F m 
itself, because sampling that kernel, put into concrete in App. El demands calculation of a 
complex-valued inverse trigonometric function, which is a costly numerical task. 
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FIG. 12: Contour levels of the number of valid decimal digits d = 2 ... 9 of F m (z), if n = 16 or 30 

sampling points are used in Salzer's integration of the inverse Laplace transform Q5U|) . 
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(n) 



I 0.00837 17061 78265 
3 0.01725 03391 11534 
5 0.02522 57393 22044 
7 0.03228 05129 34890 
9 0.03823 26639 94284 

II 0.04289 02498 71345 
13 0.04609 15264 54310 
15 0.04772 17782 62356 



71876 67547 
01977 09174 
57375 46273 
49074 48769 
05059 06029 
82958 76362 
63304 47216 
94180 43787 



1334 - 0.03493 
7117 - 0.03535 
8919 - 0.03348 
3684 - 0.02985 
0158 - 0.02477 
4358 - 0.01853 
7673 - 0.01146 
9660 - 0.00387 



88515 18794 47953 
09659 30121 29557 
90948 21254 18355 
48245 41582 95432 
02071 79697 66742 
56403 81264 01685 
24728 96512 75189 
81037 55474 09447 



06804 3103i 
12679 1042i 
16245 21511 
08487 9898i 
19368 1785i 
00724 9039i 
11368 8442i 
26721 4718i 



A 



(n 



I -(2)7.46675 12193 45759 50393 85910 48 + 
3 (3)2.91507 59384 65429 08402 81547 90 - 
5 (5)8.32343 31208 36870 55687 37470 24 + 
7 -(7)1.12187 25580 46183 78092 28439 34 - 
9 (7)5.84396 38920 01078 49663 38594 67 - 

II -(8)1.53739 97073 01945 94831 81610 90 + 
13 (8)2.10257 24343 84449 69271 33642 99 - 
15 -(8)1.04572 70576 06995 39505 45148 52 + 

TABLE IV: Extension of Salzer's 



2)2.33418 71487 56825 21567 95817 62 i 

4) 6.02533 14214 97033 76429 48817 01 i 

5) 9.23999 52597 05792 07995 48624 41 i 

6) 2.85904 20761 32552 12275 19081 14 i 

7) 1.38271 69228 73790 17106 97305 80 i 

8) 1.17150 18184 90003 20088 57634 96 i 
8)3.52009 23255 88077 30106 90523 00 i 
8)5.83415 46534 50843 20837 34946 96 i 



35j table of reciprocals of the zeros of P n (x) and Christoffel 



Numbers A^ n ' to the case n = 16. Numbers in parentheses are powers of 10 as in Tab. [IJ Half 



of the values are not shown and follow by complex conjugation: l/p\ 
(i = 2,4,6,... ,n). 



(n) 



(n) 



i-l> 
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VII. RELOCATION TO OTHER SPECIAL FUNCTIONS 



A. Gautschi's approach to the Faddeeva Function 

Following and Chapt. 7.1 of [0, F is related to the complex error function via 



F„ W = ^^>, (52) 



and to the Faddeeva function w as 



F G {z) = \J* z [l-e-*w{iVz)]. (53) 



According to 



3, (2.14)]Q(7.1.4)] (7.1.15)], w has the following representation i 



in terms 



of weights A^ and abscissae of the n-point Hermite Gauss Integration (see App. IU| : 

n > (n) 

fc=l 2 f fc 

The computation needs only about n/2 complex divisions and additions, since the weights 
and abscissa group in symmetric pairs: 

-W « T* S (55) 

Fig. 1131 considers the application with n = 20 or 32. This approach here targets the same 
region covered by Fig. El either large §lz or large $tz, but is obviously superior, since it 
first is scalable through the choice of n, and at a comparable investment into the number of 
complex operations it achieves the more accurate results. 

Eq. (2.1) in the work by Chiarella and Reichel 0| looks similar to (|55p: roughly speaking, 
the are replaced by equidistant nh, and the A^ by exp(— n 2 h 2 ). 

Strand 0| proposed a method to compute the complementary error function of z, if Qz 
is small. It starts from the presumably known complementary error function of $lz, and 
therefore is too special to be treated here. 



B. Expansion in Modified Spherical Bessel Functions 



The Confluent Hypergeometric 

rn 

Spherical Bessel Functions [ljj, (13.3.6)] 



Function in (JH) may be expanded in terms of Modified 
!Cj, which may be rewritten with 0, (10.2.24)] to 
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FIG. 13: Contour levels of the number of valid decimal digits d = 5 ... 19 of Fq(z), if the order of 
the Hermite Polynomial in (|54|) is 20 or 32. 



give 

« + l; 2*) = e^f>r(i + n) (1 " a) " 

^ n=0 ^ ' 



n=0 



/ 1 . \ -<n+l/2 

(1 + a)„ 

z z z n ( -- 
(1 + a) n \zdz 



4+1/2(2) (56) 



E( _ r(2n + 1) (i^^(I^"!^ (57) 
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n znnd\ n smhz 

V z dz ) z 



[A-l]/[2z] 

1 [z(A + 1)-(A- l)]/[2z 2 ] 

2 [(z 2 + 3)(A-l)-3z(A + l)]/[2z 3 ] 

3 [(z 3 + 15z)(A + 1) - (6z 2 + 15)(A - l)]/[2z 4 ] 

4 [(z 4 + 45z 2 + 105)(,4 - 1) - (Wz 3 + W5z)(A + l)]/[2z 5 ] 

5 [(z 5 + 105z 3 + M5z)(A + 1) - (15z 4 + 420z 2 + 945)(A - l)]/[2z 6 ] 

6 [(z 6 + 210z 4 + 4725z 2 + 10395)(^ - 1) - (21z 5 + 1260z 3 + 10395z)(A + l)]/[2z 7 ] 

TABLE V: Complexity evaluating the terms in (|57|) for small values of n, with A = exp(2z). 

with a = m + 1/2. Note that a factor (6 — a — | +n) was missing in (13.3.6) of early editions 
of Q]. The individual terms are 

271 ( \n( l , ( l - a )n T (-)" (l-0) n ^ n ^ n+1 n+2 




T<-»" U + "j ITT^^" 2 '^ K (2^ITT^ (2 " + » + 0( - »' (58) 

which indicates (i) that a direct implementation based on the formulas of Tab. IVlmav suffer 
from severe cancellation of digits if \z\ is small, and (ii) that the recurrence relations 0, 
(10.2.12)] must be used in the downward direction, for example as outlined in 0|. Table fVl 
shows that already for a small number of terms used to approximate the series, a considerable 
number of complex polynomials must be computed. Fig. El indicates that the convergence 
of the series is good close to the origin of the complex plane (explained by the fact that the 
lowest order terms of the Taylor series (1581) are 0(z n )). 

n 

The assessment follows the conclusion in |3J for a similar expansion that was a hybrid of 
the Dawson and the error function for real z: Since the evaluation of the Bessel Functions 
would be of similar complexity as a straight-forward power series for small \z\, and since the 
convergence is slow for large \z\, this type of ansatz is not competitive. 



C. Dijkstra's continued fractions 



According to Dijkstra |42j, the constraint that (0) and © serve well only inside domains 
of small or large modulus may be overcome by use of an auxiliary function K that allows a 
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FIG. 14: Contour levels of the number of valid decimals d = 3 ... 12 of F m (z), if the series (|57[) is 
truncated after the term n = 7, for m = 5 and m = 1. 



suitable continued fraction 

i-Fi(a; b + 1; z) 1 z(b + 1 — a) z(b + n — a) 



K (a, b, z) _ 

6i.Fi (a; 6; z) 

Application to (jlj) with \Fi(a; a; —z) 

1 



b + z- b+l + z- 
e~ z proposes 



b + n + z— 



(59) 



F m{z) = -e z K(a, a, -z), (a = m + 1/2). 



(60) 
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In Fig. we investigate the readiness of the representation. It is excellent nearby the 
negative real axis, staying above d = 10 (m = 3,N = 16), d = 13.1 (m = 1,N = 32) and 
d = 16.6 (m = 3, N = 32). Dijkstra 42j demonstrates that K(a,b,z) mediates between 
low- z and high-z expansions for real positive z. In so far, the sign change of z in (J6U|) means 
our plot actually looks at the "wrong" side. The continued fraction is terminated at the 
iVth convergent. About one decimal in accuracy has actually be gained by deletion of the z 
in this closing denominator like 

. * 1 z Nz , . 

K(a, a,z) = — . (61) 

v / a + z- a + l + z- a + N v/ 



A combination of 



19t (13.4.4)], the Kummer transformation (13.1.27)] and the defi- 



nition (JoU)) enforces that both arguments of z carry the same sign: 

zK(l, a,z) = l- e ' (62) 

(2m - l)F m -i{z) 

Fig. El shows that this indeed extends the fitness of this Dijkstra representation into the 
region of 5fc > 0, at the expense of the fitness to compute F m (z) at $lz < 0. 
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FIG. 15: Contour levels of the number of valid decimal digits d = 1 ... 16 of F m (z), if Dijkstra's 

representation is used with (|fil|) truncated at a depths N . 
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FIG. 16: Contour levels of the number of valid decimal digits d = 1 ... 18 of F m (z), if Eq. (|62|) is 

used with Dijkstra's continued fraction (|59jl truncated at a depths N. 
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VIII. SUMMARY 



Working with standard numerical methods on the fundamental integral representation of 
the Incomplete Gamma Function F m (z) is generally inefficient as it demands dense sampling 
(frequent evaluation) of complex exponentials. 

Continued fraction and rational function approximations are difficult to control, because 
the regions of known accuracy in the z-plane are complicated. 

Being an analytic function of z, the fastest evaluation uses Taylor series which recall 
tables of the derivatives d n F m (zo) / dz$ that have been computed off-line to high precision 
(Sec. EE)). 
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APPENDIX A: EXAMPLE APPLICATION: TRAVELLING ORBITAL S 



Let \l/o( r , t) = exp(iu!t)(po(r) be a solution of the Schrodinger equation 



-^V 2 r + \/(r) 



d 

#o(r,t) = -ih—* (r,t) 



(Al) 



associated with a stationary potential V(r). The transition to a potential that moves with 
constant velocity v in the laboratory system is useful to describe electrons bound to scat- 
tering atoms, and reads 



h 2 

V 2 r + V(r - vt) 

2m 



d 

* v (r,t) = -ih— tf v (r,t) 



The solution 



h 2 k 2 

^vfr, t) = expUiut) exp(i 1) exp(— ik ■ r){p Q (r — vt) 

2m 



(A2) 



(A3) 



with hk = mv is generated from the solution of (|A1|) [43]. This Galilean transformation 
lets "travelling" orbitals oc exp(— ik ■ (r — R))y n(r — R) become a natural choice for basis 



functions in the laboratory coordinate system [44|. If the (po are linear combinations of 
Gaussian Type Orbitals (GTO's), the Coulomb Integrals may be treated with the product 



rule for GTO's 



ABO 



15 



45 



46 



47 



48 



49 



501 ]. and the Gauss transform 



1 

r 



e~ r2s2 ds 



(A4) 



of the Coulomb potential eventually reduces the Coulomb integrals to Incomplete Gamma 
Functions with complex argument z j^l], 0, 13] . 



APPENDIX B: LAPLACE REPRESENTATION 

Laplace transformation of the power series ^ on a term-by-term basis yields 
(7.621.4)] 



1 1 



(-) 



ra+1 



u n=l v ' 

oo 



;U-2m 



1 



\2m+2n-l 



(^i)i-2m ^ 2m + 2n - 1 

\ V / 71=1 V V ' 

■r- (m »^^)^ (a ';i |iH/ 



(Bl) 
(B2) 
(B3) 
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where a = m + 1/2. This result could also be derived from the Laplace transform of the 
differential equation — which is a combination of (j4T)j) and (jUJ) — : 



d . . 



-zF m+ i(2;) 



-aF m (z) + — , 



- A [ a (£F m )(a) - F(z = 0)] = -a(£F m )( S ) + 



This inhomogeneous differential equation 



S ^(£F m )( S ) = (a - l)(£F m )( S ) - l-L. 



(B4) 
(B5) 

(B6) 



is solved by writing down the solution of the separable homogeneous differential equation, 
CF m (s) = s a ~ 1 c, then introducing the function c(s) for the constant c, which leaves a simple 
differential equation for dc(s)/ds and 

1 



c(s) 



s a (s + 1) 



-ds. 



(B7) 



This is solved with the recursion j3; (2.249)], because a is half integer, until j3 5 2.211] is 
applicable for the reduced remnant. 

To round off this excursion: One could introduce (4.4.42)] 

J2 \ 2 



arctan a 



1 + a 2 



2 a 2 

1 + 



2-4 



31 + a 2 3-5Vl + a 



H h 



(2n)!! 



(2n 



a: 



+ 



(B8) 



in (jB3|) . the simplest case of m = reading 

1 



£(*b(s)) 



1 + s 



1+?- s 



31 + s 



+ 



(2n)!! 
(2n+ 1)!! 



+ 



7T 1 . . 

+ 27S- (B9) 



Truncation of this series after some nth term followed by an inverse Laplace transform 
(29.2.21)] yields e~ z multiplied by a polynomial of degree n in z, plus a/vt/z/2, i.e., replaces 
w(iy/z) I ' yfz in (jHHjl by a polynomial in z. Obviously, that is a poor representation, as it 
forces the polynomial in z to compensate the singularity at z — 0. 



APPENDIX C: ROOTS OF HERMITE POLYNOMIALS 

These pairs of weights and abscissa of the Gauss-Hermite quadrature can be taken from 
[19I (Tab. 25.10)] for some values of n < 20, for n = 8, 16, 32 or 64 from j^J, or otherwise 
computed with the dOlbcf routine of the NAG library. Following on a reference by Shao 
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et al. 3| to a note by Hofso miner |55j |. the zeros of H n (x) can be refined through a third 
order Newton method. Improved solutions x^ +1 > are computed from guesses x^ through 



5jJ (4)], here 



n 



H n (x^) 

HL(xM) 



x 



H^(xO))\ • 



(CI) 



Note two sign errors in 55, (3)]; the correct equation is 

« = x-f/f'-(P + S/f')(f/f) 2 
1 



(4 p; - P> + Q + lopfif//' - S"//' + 6S 2 /f 2 )(f/fy + 0[(/// 



p/\4l 



(C2) 



Similar to [5J, (5.10)] one might consider using the terminating continued fraction 
H' n {x) 2n 2(n - 1) 2{n - 2) 2 n (n - 1) (n - 2) (n - 3) 



H„ (x) 2x- 



2x- 



2x- 



2x x— 2x- 



x— 2x- 



(C3) 



in (jClj) to meet the thread of cancellation of digits. The weights follow as (25.4.46)] 



A 



(n) ^ 2 w - 1 ra!0F 



(C4) 



APPENDIX D: "PERTURBED" QUADRATURE 







Eq. (J47|) is an application of the Euler-Mclaurin formula |l2l (25.4.7)] 



(Dl) 



This is in contrast to a rule that involves the value in the midpoint of the integration intervals 
d (2)]: 



f(x)dx =(b- a)f [ Z±*) + [f{b) - f{a)\ . (D2) 



2 ) 24 

The evaluation of the function at the end points a and h in Sect. IVI El is cheaper than an 
additional evaluation in the middle of the integration interval, because it needs two complex 
multiplications but no new exponentials. So the midpoint rule and variants proposed by 
Hammer and Wicke (5?], 0, 59, Q are not advantageous in our case. A cubic spline 
interpolation f(x) = Yln=o Cn%n induces the higher moments 

b — a 



xf(x)dx 



60 



a 2 [2f(b) - 3f'(a)] + a [&/'(&) + 9/(6) + 21/ (a) + bf(a)} 



2b 2 f(a) + 2lbf(b)+9bf(a)-3b 2 f(b) 



(D3) 
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x 2 f(x)dx = ^Jt [a 3 if{b) - 2f(a)} + a 2 [6f (6) + 4/(6) + 16/(o)] 
+ab [bf'(a) + 10/(6) + 10/(o)] 

+b 2 [4/(a) - 2b f '(b) + 16/(6) + bf'(a)} }, (D4) 

which establishes (14^1). and 



x 3 f(x)dx = b —£ {a 4 [4f (6) - 10/'(a)] + a 3 [56/'(6) + 90/(a) - 2b f (a) + 15/(6)] 



420 



+3a 2 6 [&/'(&) + 13/(6) + bf'(a) + 22/(a)] 
+a6 2 [66/(6) + 39/(a) + 56/'(a) - 26/'(6)] 

+6 3 [-106/'(6) + 46/'(a) + 15/(a) + 90/(6)] }. (D5) 



V/(x)dx = ^ {a 5 [5f (6) - 15/'(o)] + a 4 [76f (6) + 150/(o) - 5b f (a) + 18/(6)] 

+2a 3 6 [36/'(6) + 24/(6) + bf'(a) + 60/(a)] 
+2a 2 6 2 [42/(6) + 42/(a) + 36/'(a) + bf{b)\ 
+ab 3 [-56/'(6) + 48/(a) + 76/'(a) + 120/(6)] 
+6 4 [-156/' (6) + 56/'(a) + 18/(a) + 150/(6)] }. (D6) 

The generalization of ()D1|) to an integral over a quintic spline that engages also the second 
derivatives (curvatures) at the interval limits reads 

b f(x)dx =(b-a) M±M _ [f{b) _ f{a)] + (^1! [f{b) + f (a) ] ; (D7) 

which simplifies to (|D1|) if f(x) is any cubic polynomial. The three formulas ()D5j) - (jD7j) have 
not been used in this work. 



APPENDIX E: GAUSS-JACOBI ABSCISSAE AND WEIGHTS 

Formalas of weights Wi and abscissae t« with (jUJ) are given in Q, (25.4.33)], where ti are 
the zeros of Jacobi Polynomials Pn 2m '°\l — 2t), and 

n-1 

1/wi = J>m + 2 J + l)[^ (2m ' 0) (l - 2t,)] 2 - (El) 
i=o 

I m — 0, this reduces to the Gauss-Legendre quadrature, Table 25.4 in |l9j. Table 25.8 in 
I9I ] covers the cases m < 2 with n < 8, and we provide Tables IVHIVHI to cover 2m = k = 2 
or 4 with n = 20. 
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W; 



1 


(-1)0.14204 21115 93581 53319 97686 86 


(-5)0.37492 


20993 33713 47725 


24136 


88 


2 


(-1)0.37851 28784 95018 10290 83919 36 


(-4)0.41039 


10208 73202 05549 


07404 


17 


3 


(-1)0.71300 98508 


12494 54654 90216 28 


(-3)0.19388 


83096 17511 81076 


98385 


25 


4 


(0)0.11385 


86970 


85452 


85631 


27170 


70 


(-3)0.60704 


57042 67258 25415 


27237 


40 


5 


(0)0.16462 


10853 


68438 


27482 


84816 


05 


(-2)0.14774 


44441 97148 58062 


83857 


64 


6 


(0)0.22250 


74079 


64513 


24644 


88740 


54 


(-2)0.30224 


89170 90737 93749 


18485 


98 


7 


(0)0.28628 


43889 


84712 


31925 


79165 


65 


(-2)0.54320 


90806 12302 74500 


12011 


14 


8 


(0)0.35459 


29542 


22632 


43338 


26791 


65 


(-2)0.88135 


56078 24335 53289 


25273 


17 


9 


(0)0.42597 


73418 


17173 


12827 


24878 


03 


(-1)0.13140 


91925 66616 51742 


10813 


37 


10 


(0)0.49891 


61845 


17151 


44862 


93753 


46 


(-1)0.18220 


50341 82686 91551 


84331 


53 


11 


(0)0.57185 


49590 


66743 


31538 


45749 


65 


(-1)0.23682 


29041 42460 53980 


86645 


03 


12 


(0)0.64323 


91298 


97546 


89435 


29449 


33 


(-1)0.29002 


40474 55892 85169 


91635 


73 


13 


{{J JU. ( 11D4 


( loio 


1 Q707 






yo 


(-1)0.33556 43145 34754 9347900271 


48 


14 


(0)0.77532 


35806 


15694 


96159 


06180 


35 


(-1)0.36697 42047 20181 32713 25424 


96 


15 


(0)0.83320 


87465 


18990 


56998 


70192 


99 


(-1)0.37847 39037 69223 40715 35469 


55 


1 

1U 


(0)0.88396 


90923 


44513 


19505 


33025 


68 


(-1)0.36587 91655 32270 64986 60929 


so 


17 


(0)0.92652 


28082 


14714 


71694 


24334 


55 


(-1)0.32734 64431 75731 48047 47964 


86 


18 


(0)0.95996 


30995 


38093 


21900 


36798 


50 


(-1)0.26382 56442 91255 84794 51338 


97 


19 


(0)0.98357 


79118 


66012 


16709 


74583 


86 


(-1)0.17913 75232 57385 59101 77395 


77 


20 


(0)0.99686 


93162 


59256 


41043 


78498 


58 


(-2)0.79757 92736 27665 16852 27309 


03 



TABLE VI: Extension of [19|, Tab. 25.8] and |6J| to n = 20, k = 2: the abscissas and weights for 



the Gaussian integration of moments, L x k f(x)dx « 2^=1 w if( 
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X j 



(-1)0.28236 72218 29331 58389 89347 65 (-7 

(-1)0.59393 81548 17751 78032 25227 62 (-6 

(-1)0.98825 08311 50225 26221 23471 31 (-5 

(0)0.14592 04712 27084 45276 93121 00 (-4 

(0)0.19980 20982 26471 39701 42744 40 (-4 

(0)0.25943 51729 24709 24862 43232 70 (-3 

(0)0.32366 51708 61772 75707 78323 99 (-3 

(0)0.39124 50422 79918 49074 92563 95 (-2 

(0)0.46086 11522 78328 67423 18689 83 (-2 

(0)0.53115 95358 51314 54342 32325 30 (-2 

(0)0.60077 25678 99009 31814 54529 39 (-2 

(0)0.66834 57418 12982 73134 52276 12 (-1 

(0)0.73256 41195 17517 70122 33303 84 (-1 

(0)0.79217 79750 00657 56854 95189 99 (-1 

(0)0.84602 71502 09447 61169 84320 65 (-1 

(0)0.89306 36602 60154 22796 45053 63 (-1 

(0)0.93237 21209 85837 16946 52904 79 (-1 

(0)0.96318 76416 89199 08774 24334 04 (-1 

(0)0.98491 10827 62489 56330 49229 13 (-1 

(0)0.99712 45845 24283 68493 64961 69 (-2 

TABLE VII: Extension of |ig|, Tab. 25.8] 



)0.17148 95670 15906 66781 55999 02 
)0.44002 56392 91105 57895 18288 23 
)0.41390 07241 19088 67017 87192 60 
)0.22963 59798 81845 93019 79823 98 
)0. 90745 14619 04243 16233 00990 36 
)0. 28147 15320 39879 02977 22252 87 
)0.72562 36166 25890 38003 96362 88 
)0.16125 65884 61398 80639 52699 59 
)0.31660 97873 98703 11170 19225 67 
)0.55864 33047 22474 38067 86431 36 
)0.89646 32789 96732 44746 29624 35 
)0.13190 88993 08882 52772 58617 37 
)0. 17889 49674 56314 16340 77879 83 
)0.22414 41002 15713 19313 92980 20 
)0.25926 78294 51372 33371 02635 30 
)0.27551 70856 12976 95932 96393 76 
)0.26583 07034 97369 16858 79230 15 
)0.22683 18042 40368 25551 58049 92 
)0.16017 53993 12516 57464 20909 00 
)0.72877 91419 97403 30297 32866 81 

and Tab. EUto n = 20, k = 4. 
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